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ABSTRACT 

Self- Interacting Dark Matter (SIDM) is a collisional form of cold dark matter (CDM), 
originally proposed to solve problems that arose when the collisionless CDM theory 
of structure formation was compared with observations of galaxies on small scales. 
The quantitative impact of the proposed elastic collisions on structure formation has 
been estimated previously by Monte Carlo TV-body simulations and by a conducting 
fluid model, with apparently diverging results. To improve this situation, we make 
direct comparisons between new Monte Carlo iV-body simulations and solutions of 
the conducting fluid model, for isolated SIDM haloes of fixed mass. This allows us 
to separate cleanly the effects of gravothermal relaxation from those of continuous 
mass accretion in an expanding background universe. When these two methods were 
previously applied to halo formation with cosmological boundary conditions, they 
disagreed by an order of magnitude about the size of the scattering cross section 
required to solve the so-called 'cusp-core problem.' We show here, however, that the 
methods agree with each other within 20 per cent for isolated haloes. This suggests 
that the two methods are consistent, and that their disagreement for cosmological 
haloes is not caused by a breakdown of their validity. 

The isolated haloes studied here undergo gravothermal collapse. We compare the 
solutions calculated by these two methods for gravothermal collapse starting from 
several initial conditions, including the self-similar solution by Balberg, Shapiro, & 
Inagaki (2002, BSI), and the Plummer, Navarro-Frenk- White and Hernquist profiles. 
We compare for the case in which the collisional mean free path is comparable to, or 
greater than, the size of the halo core. This allows us to calibrate the heat conduction 
which accounts for the effect of elastic hard-sphere scattering in the fluid model. The 
amount of tuning of the thermal conductivity parameters required to bring the two 
methods into such close agreement for isolated haloes, however, is too small to explain 
the discrepancy found previously in the cosmological context. We will discuss the 
origin of that discrepancy in a separate paper. 

Key words: cosmology: theory - dark matter - galaxies: kinematics and dynamics 
- galaxies: haloes - methods: numerical - methods: iV-body simulation. 



1 INTRODUCTION 

In the currently standard cosmological model (ACDM), a 
flat universe with a cosmological constant contains colli- 
sionless cold dark matter as its dominant matter compo- 
nent, perturbed by primordial Gaussian-random-noise den- 
sity fluctuations. This model has been highly successful 
at explaining observations of the background universe and 
large-scale structure. On small scales, however, the distribu- 
tion of dark matter in coordinate and velocity space is not 
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fully understood. A-body simulations show that the den- 
sity profiles of the virialized regions ('haloes') that form in 
this collisionless dark matter are cuspy, such as the Navarro, 
Frenk, & White (1997, NFW) profile, in which p -> r _1 to- 
ward the centre. Recent high-resolution simulations show 
that the inner profile is not exactly a power law (Hayashi 
et al. 2004; Diemand et al. 2005; Merritt et al. 2006; Gao 
et al. 2008; Stadel et al. 2009; Navarro et al. 2010), but 
it diverges nevertheless. On the other hand, a cored profile 
(a profile which flattens in the centre) such as a pseudo- 
isothermal profile, is favored by observations of dwarf and 
low surface brightness (LSB) galaxies (de Blok et al. 2003; 
Gentile et al. 2004; Kuzio de Naray et al. 2006; Zackrisson 
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et al. 2006; Gentile et al. 2007; Kuzio de Naray et al. 2008; de 
Blok et al. 2008; Oh et al. 2010b; de Blok 2010). Dwarf spi- 
ral and LSB galaxies are dark-matter dominated. As such, 
it was originally thought, their mass distribution should re- 
flect the dark matter dynamics alone, and be relatively less 
affected by the complexity of the dissipative baryonic com- 
ponent. This, it was thought, makes these systems ideal for 
studying the undisturbed, intrinsic, dark matter distribution 
on small scales. 

This apparent cusp-core conflict was one of the small- 
scale structure problems of the CDM model which prompted 
suggestions that the dark matter might be something else, 
with microscopic properties that would alter the structure 
on small scales without spoiling the success of CDM on large- 
scales. Spergel & Steinhardt (2000) proposed self-interacting 
dark matter (SIDM) as a possible solution to this cusp-core 
problem, adding hypothetical elastic-scattering collisions to 
the otherwise collisionless particles of the standard CDM 
cosmology. Heat transfer within the virialized haloes, due 
to these non-gravitational collisions, was then suggested to 
make the halo cores expand. This latter idea was confirmed 
by several numerical and analytical studies. Burkert (2000) 
introduced a Monte Carlo scattering algorithm between dark 
matter particles to take the self interaction into account 
in a numerical JV-body simulation, and this method was 
refined by Kochanek & White (2000). These iV-body re- 
sults suggested, however, that SIDM haloes would undergo 
gravothermal collapse, making them unsuitable to explain 
the observed haloes. 

Balberg, Shapiro, & Inagaki (2002, hereafter BSI) ap- 
plied a conducting fluid model, originally invented to de- 
scribe gravitational scattering in star clusters, to isolated 1 
SIDM haloes. BSI derived a self-similar gravothermal col- 
lapse solution at large Knudsen number (Kn 3> 1, where 
Kn is the ratio of SIDM scattering mean free path to the 
system size). The isolated halo collapsed within a finite time. 
They also showed that the collapse is delayed compared to 
their self-similar solution when the Knudsen number is com- 
parable to or smaller than one, because the length scale of 
energy exchange is restricted by the mean free path. Their 
conclusion that the collapse time would naturally exceed a 
Hubble time was more optimistic about the SIDM hypoth- 
esis. 

A realistic halo is not isolated, however, since it forms in 
a cosmologically expanding background universe, with infall 
and a finite pressure at the virial radius (Shapiro et al. 1999). 
Cosmological simulations showed that a cross section per 
unit mass a — 0.5 — 5 cm 2 g _1 makes the profile cored, 
and the cored profile is stable (Yoshida et al. 2000; Dave 
et al. 2001). Colin et al. (2002) emphasized that the profile 
depends on the accretion history, especially when the last 
major merger occurred. 

With the importance of cosmological infall in mind, Ahn 
& Shapiro (2005) derived a cosmological similarity solution 
for this problem, with proper account taken of such cosmo- 



1 The term 'isolated halo' shall refer here to an object of fixed 
mass with vacuum boundary conditions; i.e., not subject to cos- 
mological boundary conditions involving mass infall or evolution 
by perturbation growth in a Fricdmann-Robcrtson- Walker uni- 
verse. 



logical boundary conditions. This solution shows that the 
gravothermal collapse, which occurs for the isolated halo, is 
prevented by infall in the cosmological case, and the core 
has a constant size in units of the virial radius, for a given 
SIDM cross section. When there is no self-interaction, this 
fluid approximation gives a density profile similar to the 
cuspy profile found in iV-body simulations. This shows that 
the fluid approximation also describes the virialization of 
collisionless dark matter appropriately, providing in effect 
an analytical derivation of the NFW profile in the colli- 
sionless limit. This is because the collisionless Boltzmann 
equation reduces to fluid conservation equations for an ideal 
gas with ratio of specific heats equal to 5/3 if the velocity 
distribution of the particles in the frame of bulk motion is 
isotropic and skewless (see Section 2.2). In the presence of 
SIDM scattering collisions, too, those analytical solutions 
are in qualitative agreement with the corresponding Monte 
Carlo iV-body simulations; SIDM haloes have cores and the 
cores collapse within a finite time when they are isolated, 
but they do not collapse within a Hubble time in a cosmo- 
logical environment. However, the values of the cross section 
necessary to explain the observed dark matter density pro- 
files are not in agreement. Ahn & Shapiro (2005) find that 
a ~ 200 cm 2 g _1 fits the dwarf and LSB galaxy rotation 
curves best, while iV-body simulations suggest the range of 
a = 0.5 — 5 cm 2 g _1 gives the observed central density. 

In order to test the SIDM hypothesis by compari- 
son with astronomical observations, it is necessary to im- 
prove our understanding of these theoretical predictions 
and reconcile them. That will be the focus of this pa- 
per, as described below. In the meantime, related progress 
continues to be made on other fronts, in comparing both 
SIDM and collisionless CDM with observation. For exam- 
ple, non-circular motions may affect the density profile es- 
timate (Hayashi & Navarro 2006), but they are usually not 
strong enough to make observations consistent with the 
theoretically-predicted cuspy profile (Oh et al. 2008; Tra- 
chternach et al. 2008; Kuzio de Naray et al. 2009). The 
cusp-core conflict may be mitigated by baryonic processes 
of gas outflow, induced by supernovae feedback (Governato 
et al. 2010; Oh et al. 2010a). Central dark-matter densities 
could also be reduced, in addition, by bars (Weinberg & 
Katz 2002), or by gas clumps sinking via dynamical friction 
(El-Zant et al. 2001), while the system is baryon rich, but 
these two scenarios may not be strong enough to make large 
enough cores in realistic situations (Sellwood 2008; Jardel & 
Sellwood 2009) . While observations of dwarf and LSB galax- 
ies generally prefer cored density profiles, some, however, are 
also consistent with an NFW profile (Hayashi et al. 2004; Si- 
mon et al. 2005; Hayashi & Navarro 2006; Valenzuela et al. 
2007). The wide diversity of dwarf/LSB cores has also led 
to a suggestion that SIDM alone cannot be the full solution 
to the cusp-core problem (Sanchez-Salcedo 2005; Kuzio de 
Naray et al. 2010). 

There are also some constraints on the value of the 
SIDM cross section from comparison of theoretical predic- 
tions with galaxy clusters, for which dark matter velocity 
is much higher than for galaxies. A r -body results find that 
the core size of relaxed SIDM clusters of galaxies becomes 
too large for cross section a > 1 cm 2 g _1 (Yoshida et al. 
2000; Arabadjis et al. 2002; Lewis et al. 2003). The ana- 
lytical cosmological self-similar solutions of Ahn & Shapiro 
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(2005), however, point out that the core sizes of clusters are 
also small enough, not only for small cross section in the 
long-mean-free-path regime, but also for large cross section 
in the short-mean-free-math regime, i.e., a > 200 cm 2 g _1 , 
because the short mean free path then limits the amount of 
heat conduction. However, such a large cross section would 
also enhance the fluid-like behaviour of SIDM, which may 
then conflict with observations of merging clusters. For ex- 
ample, observations of cluster IE 0657-56, known as the 
'bullet cluster,' from which total matter density has been 
mapped by weak and strong gravitational lensing measure- 
ments while the density of the intergalactic baryon-electron 
fluid was mapped by X-ray measurements, show that the 
dark matter spatially segregated from the baryon-electron 
plasma, as it would be if the dark matter is not highly colli- 
sional (Clowe et al. 2006). The dark matter and galaxies of 
the subcluster have passed through the main cluster with- 
out distortion while the baryon gas shows a bow shock due 
to its collisional nature. This observation excludes the pos- 
sibility that SIDM is too highly collisional, or fluid-like. An- 
alytical estimates and Monte Carlo iV-body simulations of 
the bullet cluster by Markevitch et al. (2004) and Randall 
et al. (2008) constrain the velocity independent cross sec- 
tion to be a < 0.7 cm 2 g _1 , using the fact that the mass- 
to-light ratio of the subcluster is normal. Larger cross sec- 
tion would make the mass-to-light ratio smaller because the 
SIDM collisions scatter the dark matter out of the subhalo. 
[On the other hand, see Mahdavi et al. (2007) for a possi- 
bility that SIDM with cross section ~ 4 cm 2 g _1 explains a 
cluster Abell 520 which has substructures with anomalous 
mass-to-light ratio]. If the cross section is velocity depen- 
dent, it puts a constraint only at very large relative velocity. 
The relative velocity of the merging haloes is estimated to 
be between 2,500 and 4,000km s _1 at their observed sep- 
aration, and even higher at centre passage (Milosavljevic 
et al. 2007; Springel & Farrar 2007; Mastropietro & Burk- 
ert 2008). If the cross section is too large, the SIDM halo 
could also be too spherical compared to the observed ellip- 
tical matter distribution of galaxy clusters (Miralda-Escude 
2002). In addition, large cross sections result in too much 
heating and evaporation of substructure haloes in clusters 
(Gnedin & Ostriker 2001). These constraints almost ex- 
clude the possibility that a velocity- independent cross sec- 
tion solves the cusp-core problem for dwarf/LSB galaxies, 
but, velocity- dependent cross section can still be effective on 
dwarf scales, and simultaneously negligible on galaxy cluster 
scales: e.g., inversely proportional to relative speed (Firmani 
et al. 2000; D'Onghia et al. 2003), or short range Yukawa 
interaction, whose scattering cross section has v~ 4 velocity 
dependence for large relative velocities (Koda 2009; Loeb & 
Weiner 2010). 

There was an additional motivation for the SIDM hy- 
pothesis when it was first put forth, involving the overabun- 
dance of subhaloes in CDM iV-body simulation compared 
to observations of the Local Group. The number of sub- 
structures in the collisionless CDM model has previously 
been thought to be about an order of magnitude larger than 
the observed number of dwarf galaxies in the Local Group. 
Although the number of substructures is reduced by SIDM 
stripping, D'Onghia & Burkert (2003) claim that SIDM does 
not solve this problem; the cross-section that makes the pro- 
file cored (0.6cm 2 g _1 ) is apparently not efficient enough to 



reduce the number of satellites down to the observed level. 
Recent findings of nearby faint galaxies improve the agree- 
ment between observations and collisionless CDM simula- 
tions, but still require some feedback or reionization effect 
that suppresses star formation in low-mass haloes (e.g. Ko- 
posov et al. 2009; Wadepuhl & Springel 2010, and references 
therein). 

Apart from the original motivation for proposing SIDM 
as a solution of small-scale structure problem of the collision- 
less CDM model, interest in the SIDM hypothesis remains 
strong for other reasons, as well. The fundamental nature 
of dark matter is still unknown. Anomalous cosmic ray de- 
tections, reported by PAMELA, FERMI and ATIC, have 
recently stimulated new particle physics models for dark 
matter that result in the scattering interaction of SIDM as 
a secondary effect (Arkani-Hamed et al. 2009; Buckley & 
Fox 2010; Feng et al. 2010). Astronomical constraints on 
the SIDM cross section and its velocity dependence con- 
tinue to be of importance, therefore, in order to constrain 
such models. 

Further progress along these lines requires that we re- 
duce the uncertainties in the theoretical productions for 
SIDM. As described above, the quantitative estimates by 
Monte Carlo iV-body simulation and the conducting fluid 
model, of cross section values that are consistent with ob- 
servations of dwarf galaxy rotation curves, are in strong dis- 
agreement with each other. As long as this disagreement is 
unresolved, many of the additional constraints on the SIDM 
hypothesis described above, which also depend upon the va- 
lidity of the iV-body results or related analysis, will remain 
suspect. It is possible that either the iV-body or the fluid 
model does not describe the system correctly. To remove 
such possibility, we directly compare the two methods in 
the simplest case, namely, the isolated spherically symmet- 
ric halo. 

We will summarize the basis for the conducting fluid 
model in Section 2. In Section 3, we will describe our Monte 
Carlo numerical algorithm for the SIDM elastic scattering 
interaction in TV-body simulations. We test those two meth- 
ods against each other for isolated haloes in Section 4. The 
impact of our comparison results on the cosmological simi- 
larity solution by Ahn & Shapiro (2005) is discussed in Sec- 
tion 5. Finally, our results are summarized in Section 6. 



2 THE CONDUCTING FLUID MODEL 

2.1 Background: Gravothermal Collapse in Star 
Clusters 

The conducting fluid model (also known as the gaseous 
model) was first developed to study gravothermal collapse 
in globular star cluster systems, and has been shown to be 
successful in those systems. Lynden-Bell & Eggleton (1980, 
hereafter LBE) proposed a thermal conduction formula for 
collisional, gravitationally-bound systems based on dimen- 
sional analysis, and derived an analytical solution that de- 
scribes the self-similar collapse of a star cluster. That self- 
similar solution appears in the late stage of collapse in a 
Fokker-Planck calculation (Cohn 1980) and in A"-body sim- 
ulations (Makino 1996; Spurzem & Aarseth 1996; Baum- 
gardt et al. 2003; Szell et al. 2005). When the time evolu- 
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tion of a Plummer sphere is solved numerically by integrat- 
ing the partial differential equations of the fluid model, the 
resulting collapse time agrees with that from other methods 
(Goodman 1987; Heggie & Ramamani 1989). SIDM haloes 
and globular clusters are both 'collisional,' self-gravitating 
systems, but the angular distribution and velocity depen- 
dence of the collisions are different in the two cases. Stars 
obey Rutherford scattering, which is dominated by small- 
angle scattering and small velocity encounters, a oc v~ 4 , 
while the SIDM cross section we explore in this paper is 
isotropic and velocity independent. It is possible, however, 
that the SIDM interaction also obeys Rutherford scatter- 
ing via a 'dark-photon' interaction (Ackerman et al. 2009; 
Feng et al. 2009). Gravitational Rutherford scattering be- 
tween dark matter particles are negligible unless they are 
all 10 s - 1O 6 M mass black holes (Jin et al. 2005). 



2.2 The Basic Equations 

In this section, we review the conducting fluid model de- 
veloped by LBE and BSI before we compare its results 
with our simulations. As in these papers, we shall here re- 
strict the conducting model to the case of particle distri- 
butions that are spherically symmetric, isotropic in velocity 
dispersion, and quasi-static. The quasi-static approximation 
means that, while the fluid evolves thermally, it always sat- 
isfies hydrostatic equilibrium at each moment. For the prob- 
lem at hand, it is a good approximation, because the collapse 
timescale is always much longer than the dynamical time. 
The conducting fluid model is not, in general, restricted 
to quasi-static systems, however (Bettwieser 1983; Ahn & 
Shapiro 2005). Deviation from isotropic velocity dispersion 
is also possible to consider, but it plays only a minor role 
for the problem of interest here. 

The Boltzmann equation, which is a partial differential 
equation in phase space, can be written as an infinite series 
of moment equations in position space by integrating over 
all velocities. When third moments are negligible (skewless), 
the series of moment equations can be closed by truncating 
at second order, as described in Ahn & Shapiro (2005). The 
equations which result in this case are what we shall here 
refer to as 'the fluid approximation.' The same equations 
would result from assuming a Gaussian velocity distribu- 
tion (not necessarily isotropic or Maxwellian), an approxi- 
mation called 'Gaussian closure' (Levermore 1996, and ref- 
erences therein). If velocity isotropy is imposed in addition, 
the fluid approximation gives the familiar Euler equation for 
an ideal gas with the ratio of specific heats, 7 = 5/3. This 
fluid approximation can describe the structure formation of 
collisionless CDM, and accurately reproduce the CDM halo 
properties found in 3D iV-body simulations when applied to 
the problem of spherical cosmological infall (Ahn & Shapiro 
2005). 

When collisions are important (e.g. either gravitational 
scattering between stars or DM self-interaction), thermal 
conduction, which is a third order moment, should not be 
neglected. However, accurate evaluation of third moments is 
only successful in the small Knudsen number regime, Kn <C 
1. In this regime, the Navier-Stakes equation can be derived 
from the Boltzmann equation with the Fourier law of heat 
flux, 
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where p(r,t) is the mass density, v(r, t) is the one- 
dimensional velocity dispersion, a is the scattering cross sec- 
tion per unit mass, A = 1/pa is the mean free path, and b = 
25v / tt/32 « 1.38 is a constant 2 (Chapman-Enskog theory, 
e.g., Chapman & Cowling 1970; Lifshitz & Pitaevskii 1981). 
The local relaxation time is defined as t r (r,t) = l/(apav), 
where the constant a = yJlQ/n ~ 2.26 describes the colli- 
sion rate of particles that follow a Maxwellian distribution, 
defined in BSI. 

In the other limit, Kn 3> 1, Lynden-Bell &: Eggleton 
(LBE) found that an empirical thermal conduction formula 
with A replaced by the gravitational scale height (or Jeans 
length), H = \J v 2 /4irGp, explains the gravothermal catas- 
trophe of star clusters very well, i.e., 
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where C is an unknown constant of order unity. The scale 
height H characterizes the length scale that particles (or 
stars) orbit under the action of the gravitational force. 

BSI combined the two limiting forms of this thermal 
conduction into one as follows: 
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The first term inside the brackets is the LBE formula, equa- 
tion (2), which dominates in the large Knudsen number 
limit, A ^> H . In the other limit (A <C H), the second term 
inside the brackets dominates, and the conduction converges 
to equation (1), instead. BSI's formula is an empirical inter- 
polation between those two heat conductivity limits. BSI as- 
sumed C = b, but the exact value cannot be determined an- 
alytically. We determine the value of C by fitting the Monte 
Carlo TV-body data in Section 4.1.1. 

The conducting fluid model is a set of moment equations 
closed empirically by this conductive heat flux L/4nr 2 , as 
follows: 
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where M is the mass enclosed by radius r, and D/Dt is the 
Lagrangian derivative with respect to time. The first equa- 
tion describes hydrostatic equilibrium. The second equation 
is the first law of thermodynamics, in the form relating en- 
ergy transfer and entropy, l dQ = TdS = dE + pdV.' 

When A H, the fluid equations (4, 5) with the heat 
conduction equation (2) have a self-similar solution. In gen- 
eral, if a self-similar solution exists, if density is static as 
r — > oo and if the evolution timescale p c /p c is proportional 
to the relaxation time at the centre, t r , c (t) = t r (r — 0,t), 
then the central quantities evolve as, 



Pc(t)/ P c(0) = (l~t/t coU )- 2a/(3a - 2) , 



(6) 



2 BSI used a value b 
257I-/32. 



25^/(32^) ~ 1.002, but it should be 
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vl{t)/vl{0) = (1 - *A coii r (2Q - 2)/(3< *- 2) , (7) 

tr,c(*)/tr,c(0) = l-*/tco!J, (8) 

for some constants a and t co i;, by dimensional analysis 
(LBE). The exponents of 1 — t/t co ll depend on the form of 
the relaxation time, therefore, different from those for star 
clusters. BSI obtained, 

a = 2.190, (9) 

t coU = 290C , - 1 t r , c (0), (10) 

by solving an eigenvalue problem of a system of ordinary 
differential equations. 

For the non-self-similar time evolution considered in 
Sections 4.1.2, 4.1.3, and 4.2, we must integrate the time 
evolution of fluid variables p and v 2 numerically by alterna- 
tive steps of the heat conduction and adiabatic relaxation 
to hydrostatic equilibrium, as described in BSI. 



3 iV-BODY METHOD WITH MONTE CARLO 
SCATTERING 

3.1 Scattering Algorithm 

In this section, we will describe the Monte Carlo algorithm 
we implemented to model non-gravitational scattering of 
dark matter particles by other dark matter particles, within 
a pre-existing gravitational JV-body method. Our scattering 
algorithm is similar to Kochanek & White (2000). Each par- 
ticle can collide with one of its k nearest neighbors with a 
probability consistent with a given scattering cross section. 
For simplicity, we assume collisions are elastic, velocity in- 
dependent, and isotropic in the centre of mass frame, but 
the Monte Carlo method can handle any differential cross 
section. We first outline the Monte Carlo iV-body method 
and then explain, in detail, the algorithm that we have 
implemented in the parallel iV-body code GADGET 1.1 
(Springel, Yoshida, & White 2001), which uses the tree al- 
gorithm to calculate gravitational forces. 

Monte Carlo algorithms for particle-particle scattering 
(known as direct simulation Monte Carlo) have been used 
for more than thirty years to solve physics and engineering 
problems of collisional molecules, giving reasonable results 
(Bird 1994). For example, the results agree with an exact 
solution of the spatially homogeneous Boltzmann equation 
that describes the relaxation toward a Maxwellian distribu- 
tion; they also agree with the Navier-Stokes equation solu- 
tions and experiments, including the thermal conductivity, 
in the small Kn regime (e.g., Nanbu 1984; Bird 1994; Gallis 
et al. 2004). 

Consider iV-body particles at positions Xj and veloci- 
ties Vj with equal mass m. We discretize the distribution 
function / with, 

/(x, v) = ™W(x - x,-;rf V(v - v,-), (11) 

3 

where W(x;rk) is a spline kernel function of size ru, r*' is 
the distance from particle j to its k « 32nd nearest neigh- 
bor, and 8 is the Dirac delta function. Our choice of kernel 
is often used in smoothed particle hydrodynamics, including 
GADGET. Our algorithm is identical to that of Kochanek 
& White if a top-hat kernel is used for W instead of a spline. 



The result does not depend on the details of the kernel, how- 
ever. We tested with k = 128 but did not see any difference. 

The collision rate F for a particle at position x with 
velocity v to collide with this distribution / is, 

r = ^mIF(x-x J ;rf>|v-v J |, (12) 

3 

where a is the scattering cross section per unit mass. There- 
fore the probability that an A"-body particle collides with 
particle j during a small timestep At is, 

P 0j = mW(x -x j ;r$ th )a\v -v j \At (13) 

One can generate a random number and decide whether 
this collision happens and reorient velocities when they col- 
lide. This method is similar to a variant of direct simula- 
tion Monte Carlo called Nanbu's method (Nanbu 1980). His 
Monte Carlo algorithm, with the pairwise collision probabil- 
ity, equation (13), can be derived from the Boltzmann equa- 
tion as described in his paper. Conversely, results of Nanbu's 
numerical method converge, mathematically, to the solution 
of the Boltzmann equation as the number of particles goes 
to infinity (Babovsky & Illner 1989). In Nanbu's method, 
only one particle is scattered per collision (only particle 
but not j). The philosophy is that the A"-body particles are 
samples chosen from real sets of microscopic particles, and 
those samples should collide with a smooth underlying dis- 
tribution function, not necessarily with another sampled N- 
body particle. However, then the energy and momentum are 
not conserved per collision. Moreover, the expectation value 
of the energy decreases systematically (Greengard & Reyna 
1991). In our case, the error in the energy rises by 10 per cent 
quickly, so we decided to scatter A"-body particles in pairs, 
not using Nanbu's method. Scattering in pairs is common in 
direct simulation Monte Carlo (e.g. Bird's method). 

When particles are scattered in pairs, other particles j 
can scatter particle during their timestep as well, but the 
scattering probability Poj, in equation (13), is similar to, but 
not exactly equal to Pjo, due to the difference in kernel sizes. 
Therefore, we symmetrize the scattering probability by tak- 
ing the average scattering rate. Note that it is not trivial to 
generalize the pairwise scattering algorithm to simulations 
with unequal A"-body particle masses, because Poj and Pjo 
would then differ by a factor of their mass ratio; there is no 
reason to symmetrize two intrinsically different probabilities 
into single pairwise scattering probability. 

In the following, we describe our algorithm in detail. 
Each particle, say particle 0, goes through the following 
steps, (i) to (iii), during its timestep Ato. Let particles 
1, . . . ,k be the k nearest neighbors of particle (k — 32 ±2). 
The particle collides with its neighbors with probabili- 
ties Pjo/2 (equation 13) during its timestep. The factor of 
two is the symmetrization factor that corrects the double 
counting of pairs. A particle j would also scatter particle 
during its timestep, which results in a symmetrized scat- 
tering rate. Imagine a probability space [0, 1], with disjoint 
subsets Ij = [Yl^Zj Plo/2, Ym=i -Ro/2) that represent scat- 
tering events between particles and j. We neglect the pos- 
sibility of multiple scattering in the given timestep. Particle 
collides with at most one of its neighbors. We restrict the 
timestep so that it is small enough that this approximation 
is good enough (see equation 15 below). We generate a uni- 
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form random number x in [0, 1], and scatter particles and 
j if x falls in a segment Ij, as described below. 

(i) In the large Kn regime, most of the particles do not 
collide with another particle in a given timestep. Therefore, 
we can reduce the computation by estimating the rough scat- 
tering probability first, and compute the accurate probabil- 
ity Poj only if necessary. First, we calculated an upper bound 
to the scattering probability, 

P = pav max Ato, (14) 

where p is the approximate density calculated from rg th via 
p — fcm/|7r(ro th ) 3 , and v max is the maximum speed of all the 
particles. If the generated random number x is larger than 
P, this means that x is not in any segment Ij , therefore, the 
particle does not collide during this timestep. 

(ii) If the possibility of collision was not rejected in step 
(i), we calculate the pairwise scattering probability Pjo and 
determine which neighbor the particle collides with. The 
index j of the collision partner is the smallest integer that 
satisfies x ^ X);=i -^j, i- e - x G ^ such j exists (otherwise, 
the particle does not collide with any neighbors). 

(iii) For particle pairs that collide, we reorient their ve- 
locities randomly, assuming an elastic scattering which is 
isotropic in the centre of mass frame. Isotropic random direc- 
tions can be generated by one square root operation, without 
using trigonometric functions (e.g., Vesely 2001). The veloc- 
ities are updated in the kick phase of the leap-frog time inte- 
gration in the GADGET 1.1 gravitational iV-body method. 
At that time, we also update the centre-of-mass velocities 
of the nodes around the scattered particles in the oct-tree, 
used for the gravity calculation. This is because the centre- 
of-mass velocities can be changed drastically by scattering. 

We allowed at most one scattering per timestep per par- 
ticle. In order to suppress the error due to possible multiple 
scattering, we restrict the timestep so that, 

P<0.1. (15) 

This restriction makes the Monte Carlo method computa- 
tionally costly in the small Kn regime, because then the 
timestep becomes much smaller than the dynamical time - 
of the order of the timestep in collisionless iV-body simu- 
lations. This timestep problem may seem to be avoided by 
performing multiple scatterings per dynamical timestep, but 
there is another limit when the Knudsen number is small. 
The distances to kth neighbors must be smaller than the 
mean free path A = l/(pa), 

r k < A. (16) 

Otherwise, particles more than a mean free path away would 
be allowed to collide and, thereby, make the heat transfer 
larger than it should be in the diffusion limit. If one tries to 
avoid this by choosing a kernel size smaller than the mean 
particle separation length, then the particles simply freely 
stream beyond the mean free path, which is again incorrect. 
The only way to overcome this problem is to increase the 
number of the particles in inverse proportion to the volume 
within the mean free path, -/V par ticics oc A~ 3 = (per) 3 , which 
increases very rapidly during core collapse as p increases. 
Conditions (15) and (16) prevent us from running simula- 
tions with small mean free path. 



3.2 Dimensions and Units 

We use the following characteristic scales in the rest of 
this paper. For the initial BSI self-similar profile (Sections 
4.1.1, 4.2) and the Plummer model (Section 4.1.2), we de- 
scribe densities and velocities in units of central quantities, 
p c (t) = p(0,t) and v 2 (t) = v 2 (0,t); p(r,t) and v(r, t) de- 
note the density and the one dimensional velocity disper- 
sion, respectively, in spherical symmetry. We use core radius 
r c (t) = v c / \/ '4nGp c as a standard length scale. 

For Hernquist and NFW profiles (Section 4.1.3), which 
have singularities at the centre, we use the scale radius r s 
and density po that appear in the density profiles as units, 
instead. We use vo = r s ^AixGpo as the velocity scale, which 
is similar to the definition of the core radius above. 

The relaxation times at the centre £ r , c (£) = t r (0,t) for 
BSI and Plummer model, or t r ,o = 1/ (apoavo) for NFW and 
Hernquist profiles, are used as unit time-scales. 

We express cross sections in a dimensionless way as 
a c = p c ur c and do = poor a , where a is the scattering cross 
section per unit mass. The inverses are the Knudsen num- 
bers, the mean free path in the unit of system size (r c or r s ). 
The evolution of the system is characterized by the Knud- 
sen number Kn only, not depending on the overall physical 
scale. 

3.3 Simulation Setup 

We generate the initial conditions for iV-body particle posi- 
tions and velocities randomly from the distribution functions 
using the rejection method (Aarseth et al. 1974). 3 The dis- 
tribution function of the BSI's self-similar solution and the 
NFW profile are calculated numerically by Eddington's for- 
mula (Binney & Tremaine 1987). The distribution functions 
of the Plummer model and the Hernquist model (Hernquist 
1990) have known analytical forms. We set the initial centre- 
of-mass velocity of the system of iV-body particles to zero 
by an overall boost. 

We truncate the initial profile at some radius 77, and 
put a simple reflecting boundary, which flips the direction 
of the radial velocity if particles are moving outward passed 
the reflecting boundary. We use 77 = 600 r c for the BSI 
self-similar profile, 77 = 58.5 r c for the Plummer model, and 
77 = 100 r s for the Hernquist and NFW profiles. The density 
at 77 is smaller than 2 x 10~ 7 p c , and the heat flux (L/4irr 2 
in Section 2) at 77, calculated from the equation for the 
conducting fluid model, is smaller than 0.02 per cent of its 
maximum value. Nevertheless, the collapse time is some- 
times sensitive to the position of this reflecting boundary. 
The collapse was slower by a factor of two when we first 
used truncation radius 77 = 58.5 r c for the BSI profile, even 
though this 77 is large enough to satisfy Antonov's criterion 
for gravothermal catastrophe (Endoh et al. 1997). We also 
tested 77 = 300 r c , for the BSI profile, and the collapse time 
changed by only 3 per cent compared to 77 = 600 r c . Hence, 
our choice of 77 = 600 r c is large enough to bring about a 
converged solution to high accuracy. 

We use two timestep criteria in addition to equation 

3 See also The Art of Computational Science vol. 11 by Hut, P. 
and Makino, J.; 

http://www.artcompsci . org/kali/vol/plummer/title .html 
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run name 


initial condition 


cross section 


N 


jy-max 


r f 


€ 


C 


BSI 


BSI self-similar 


(t c = 0.067 


4 x 64 3 


473 


600 r c 


0.1 r c 


0.75 


P 


Plummer 


(t c = 0.067 


4 x 32 3 


1455 


58.5 r c 


0.1 r c 


0.8 


H 


Hcrnquist profile 


(7 q = 0.16 


2 x 64 3 


1493 


100 r s 


0.03 r s 


0.9 


NFW 


NFW profile 


(t q = 0.088 


2 x 64 3 


794 


100 r s 


0.03 r s 


0.75 



Table 1. The long-mean-free-path regime: parameters used for each run. N is the number of W-body particles, 7V™ ax is the maximum 
number of particles inside the core, rf is the radius of the reflecting boundary, e is the initial Plummer equivalent gravitational softening 
length, which is reset sometimes overtime as central density p c (<) increases; i.e. e oc r c (t) oc v^/^/p^. The constant C is the LBE prefactor 
used for the conducting fluid model. 



(15). The timestep must also satisfy At ^ rj v v c (0)/a, and 
At ^ rja/\/Gp, where v c (0) is the initial one-dimensional 
velocity dispersion, a is the local acceleration, and p is the 
local density calculated from k — 32 nearest neighbors (see 
below equation 14). We choose the dimensionless parameters 
to be rj v = 0.02 and r)G = 0.005. With this choice, the 
initial conditions are static for several dynamical times when 
scattering is turned off. Energy conservation is satisfied to 
better than 1 per cent in all runs. 

The number of particles, gravitational softening length 
and other parameters are summarized in Table 1. For cases 
BSI and P, we reset the gravitational softening length e to 
0.1 r c (t) every time the central density increases by a fac- 
tor of 10, to avoid the numerical effect of softening on the 
density profile. 

We tested that the numerical scattering rate is correct, 
by counting the number of scatterings in the simulation of 
a non-singular isothermal sphere 4 and comparing it to the 
analytical rate. The scattering rates agree within 3 per cent 
when 64 3 particles are used and the isothermal sphere is 
truncated at 58.5 r c . The difference is due to the fluctua- 
tion in the randomly generated initial condition, not to the 
Poisson fluctuation in the number of scatterings. 



3.4 Analysis Methods 

We calculate the central quantities for each snapshot from 
A-body particles within a sphere of radius r c around the 
density- weighted centre of mass (Casertano & Hut 1985), 
assuming an isothermal sphere inside r c . Namely, we find a 
consistent r c iteratively that satisfies, 

r c = y/vl/^Gpc, (17) 

where p c is the central density estimated from M(r c ), the 
mass inside r c , 

p c ee 1.10 x M (r c )/^r 3 c , (18) 

and v c is the velocity dispersion inside r c . The value 1.10 in 
equation (18) is the ratio of the central density to the average 
density inside r c for the non-singular isothermal sphere. In 
this way, we can use as many particles as possible without 
systematic error for the calculation of the central quanti- 
ties. The velocity dispersion inside the core quickly becomes 
isothermal due to collisions. 



4 The solution of hydrostatic equilibrium (equation 4) with con- 
stant v(r) = v c and finite central density. 



We calculated the smoothed density and velocity dis- 
persion field at each point in space using an adaptive ker- 
nel, similar to that used in the well-known SPH method, 
with a Gaussian kernel whose size is rz2 - the distance to 
the 32nd nearest neighbor. These smoothed-particle density 
and velocity dispersion fields were then averaged over spher- 
ical surfaces on different radii r to give the radial profiles of 
these quantities. In practice, this amounts to summing the 
spherically-averaged Gaussian kernels evaluated at each ra- 
dius r (e.g., Reed et al. 2005). 



4 RESULTS 

We compare our Monte Carlo iV-body simulations with the 
conducting fluid model, first in the large Kn regime for var- 
ious initial conditions, and then in the transitional regime 
Kn ~ 1. 

4.1 Long Mean- Free-Path Regime 

4-1.1 Self-Similar Gravothermal Collapse Solution 

In this section, our Monte Carlo A-body simulation for large 
Kn, <t c (0) = Kn^ 1 = 0.067, is compared with the BSI self- 
similar solution. Fig. 1 shows the evolution of the density 
and velocity dispersion profiles. The right panels, plotted 
in self-similar variables, show that, when the iV-body par- 
ticles are initialized according to a given time-slice of the 
self-similar solution, they indeed evolve self-similarly in the 
Monte Carlo iV-body simulation, thereafter. 

Fig. 2 shows that our Monte Carlo iV-body simulation 
is in excellent agreement with the self-similar solution (equa- 
tions 6-10), by adjusting the value of one parameter, C, in 
the heat conduction equation (2). We determine the col- 
lapse time tcoii = 385t r , c (0) by fitting the time evolution of 
the relaxation time data (Fig. 2, left-bottom panel) by the 
linear function, equation (8). The prefactor of the thermal 
flux C = 0.75 follows from equation (10). The best-fitting 
power law index of v\ oc pi" 2 ^ a (equations 6 and 7) gives 
a value a = 2.22 (Fig. 2, right-bottom), which agrees rea- 
sonably well with the value of the self-similar solution, 2.19 
(equation 9). 

To test convergence, we simulated the self-similar solu- 
tion with three different numbers of particles, 2 x 32 3 , 4 x 32 3 , 
and 16x32 3 , which give the initial number of particles inside 
the core N C (Q) = 227, 376 and 1527, respectively. At the time 
of our resolution test, we used reflecting radius r/ = 58.5 r c . 
Numerical errors due to finite number of particles should 
only depend on the A-body particle density, independent of 
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Figure 1. Monte Carlo iV-body results for evolution from BSI 
self-similar initial conditions. Profiles of density (top) and velocity 
dispersion (bottom) for a c = 0.067 plotted in units of fixed, initial 
central values, as labelled (left), and in units of the time varying 
self-similar quantities (right) at t/t r ^ c (0) = 0,275,356,376 and 
383. The evolution is indeed self-similar. 




200 400 200 400 

f/f„(0) f/f«(0) 




200 400 1 10 100 



f/f,c(0) Pc/Pc(0) 

Figure 2. Comparison of ./V-body results and BSI similarity so- 
lution for the case shown in Fig. 1. Central quantities plotted 
as a function time. TV-body results are in good agreement with 
the self-similar solution (smooth curves, equations 6-9), with 
t C oll = 385r r ,c(0) or C = 0.75. Fluctuations are of order Pois- 
son noise; error bars represent Ap c /p c = At r ^ c /t r ,c = "2/\/N c , 
and A-ii^/v^ = l/\/N c , plotted for every 10 data points, where 
N c is the number of particles inside r c . 




200 400 600 800 200 400 600 



f/f, c (0) f/f„(0) 

Figure 3. Same as Fig. 2, except for simulations with different 
total numbers of iV-body particles. Density deviates from BSI 
solution when 7V C < 100. Collapse times are different from those in 
run BSI because the reflecting boundary is set at 58.5 r c . Density 
curves (left panel) are shifted left by 50, and up by factor of 2 for 
readability; three smooth curves are all the same BSI self-similar 
solution. 

the choice of boundary condition. As a result, it was not nec- 
essary to run additional convergence test for the case with 
r/ = 600 r c , used in our final runs. Other parameters are the 
same as run BSI. In Fig. 3, we plot the evolution of density 
and number of particles inside the core as a function of time. 
The three smooth curves in the left panel are the same ana- 
lytic self-similar solution (equation 6) with t coU = 705 t r>c (0). 
The results for the runs with N = 4 x 32 3 and N = 16 x 32 3 
agree very well, while those for N = 2 x 32 3 deviate from 
the other two runs systematically when N c < 100. Our run 
BSI contains 470 particles inside the core at t = (Table 1), 
which is better than the converged N = 4 x 32 3 run here. 
Other runs in the following sections have better resolution 
inside the core, because of the smaller fraction of particles 
for r > r c , resulting from the steeper decline in their density 
profiles. 

4.1.2 Plumraer Model 

We compare the Monte Carlo iV-body simulation with the 
conducting fluid model when the initial condition is the 
Plummer model in the large Kn regime. The Plummer 
model, a standard initial condition used to study gravother- 
mal instabilities, has a spherical mass distribution given by, 

M T 1 

P(r) _ 4^73 (l + r/a P 5 / 2 • (19) 

In this case, the characteristic scales defined in Section 3.2 
can be shown to be, v c = (2GMt / apt) 1 ^ / 12 and r c = 
a p ;/(3\/2)- We evolve this system according to the conduct- 
ing fluid model in the large Kn limit with the quasi-static 
approximation described in Section 2, equations (2), (4), and 
(5). The iV-body simulation is also performed in the large 
Kn regime, <r c (0) = 0.013. The time evolution for this quasi- 
static system is independent of the actual value of a for the 
large Kn regime, if the time for solutions with different a is 
expressed in units of the relaxation time. 

We plot the time evolution in Fig. 4. The fluid model 
agrees with the iV-body results reasonably well if the coef- 
ficient C is assigned a value of 0.80. This is in reasonably 
good agreement with the value of C = 0.75, found for the 
self-similar solution. The logarithmic slope (right-top panel) 
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Figure 4. Collapse of the Plummer model for <r c (0) = 0.013. 
( Top:) Snapshots of the Monte Carlo TV-body simulation, taken at 
t/t r>c (0) = 0.0,24.338.552.0,56.7,58.6, and 59.3. {Bottom:) Cen- 
tral density and relaxation time evolution as functions of time 
for both the TV-body results and the ID calculation of the fluid 
model with C = 0.8 (smooth curves). 



has a plateau at about —a, which is the asymptotic slope 
of the self-similar profile. This implies that the inner part 
is converging to the self-similar solution, with an asymp- 
totic logarithmic slope —a, which was well known in the 
gravothermal collapse of star clusters. 

4.1.3 NFW & Hernquist Profiles 

We also considered the case of an initial condition with 
spherical mass distribution given by the NFW profile 
(Navarro et al. 1997) to test the consistency of the Monte 
Carlo iV-body simulation and the conducting fluid model 
with each other for the typical halo profile seen in cos- 
mological TV-body simulations. The gravothermal collapse 
timescale is interesting as a way to put constraints on the 
SIDM cross section, because the result of gravothermal col- 
lapse is a divergent profile with par -2 ' 2 , which is not seen 
in observations. We note that gravothermal collapse is pre- 
vented by cosmological infall or major mergers, but when 
a halo decouples from cosmological growth, the halo could 
experience gravothermal collapse. The NFW profile is given 
by, 

p{r) = PO r/r s (llr/r s )* - (20) 

In addition, we perform a test with initial conditions 
based upon a Hernquist profile (Hernquist 1990), 

1 



P{r) = po 



r/r s (l + r/r s 



(21) 



which has the same inner profile, p oc r _1 , and is sometimes 
used to approximate the NFW profile. As we will demon- 
strate below, the gravothermal collapses for the NFW and 
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Figure 5. Gravothermal collapse of SIDM haloes starting from 
NFW and Hernquist profiles: The density and velocity dispersion 
profiles. (Left:) TV-body snapshots with NFW initial condition at 
t/t ri0 = 0, 197, 336, 462 and 515, and profiles of the fluid model 
when they have the same central density. (Right:) Same as left 
panels, except with Hernquist initial profile, plotted at t/t r fi = 
0,16,111,129 and 139. 



Hernquist profiles are significantly different, and it is there- 
fore, dangerous to use the results for the Hernquist profile 
(Kochanek & White 2000) to describe realistic haloes. As 
for the Plummer model, our conducting fluid model calcula- 
tions adopt the large Kn limit, and the TV-body simulations 
use a small, but nonzero, a (cf. Table 1). 

Figs. 5 and 6 show that the fluid model agrees reason- 
ably with TV-body results. In particular, the density profiles 
in the two cases evolve through a sequence which is close to 
the TV-body results, although the agreement is better if the 
following adjustments are made to the time axis of the fluid 
model solutions (6) . For NFW, TV-body and fluid model cen- 
tral density agree best when compared at times which differ 
by a small shift equal to 100t r ,o, or 20 per cent of the col- 
lapse time. For the Hernquist initial profile, the agreement 
is best if the time axis is scaled by a factor of order of unity 
which is equivalent to adjusting the value of C in the con- 
ductivity from C = 0.75 (found in the BSI run) to 0.9, which 
also differs by only 20 per cent. For both NFW and Hern- 
quist initial conditions, density profiles agree very well when 
the central densities are equal, but the values of the central 
densities differ by about 20 per cent at the maximum core 
expansion. 

We note that the initial NFW profile evolves very differ- 
ently from the initial Hernquist profile, for the same param- 
eters po and r s . The NFW run has a central density about 
three times smaller at maximum core expansion, and a col- 
lapse time about four times longer than for the H run. This 
is because the NFW profile has larger heat flux at r > r s 
due to its larger density there, which heats and expands the 
central mass more than does the Hernquist profile. 
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Figure 6. Central density versus time, starting from NFW (left) 
and Hcrnquist profiles (right). We use C = 0.75 for NFW and 
C = 0.9 for Hernquist profile. In the left panel, the origin for the 
fluid curve is shifted to the right along the time axis by 100 t r> o to 
show the agreement in the gravothcrmal collapse phase. Minimum 
density of the fluid is ±20 per cent different from .W-body. Error 
bars show Ap c /p c = 2j \[W C for every 10 data points. 
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Figure 7. Velocity anisotropy for JV-body results. The radial 
and tangential velocity dispersion for each run. The snapshots 
are taken when the densities at the centre are 10p c (0) for run 
BSI and P, 25pn for run NFW and 125pn for run H. The velocity 
dispersion is calculated in 40 equally-spaced logarithmic bins. The 
error bars show the Poisson fluctuations when they are larger than 
1 per cent, Av 2 = v 2 /\/N, where N is the number of particles in 
each bin. 



Our simulation results are qualitatively similar to those 
of Kochanek & White (2000), but the time evolutions dif- 
fer by a factor of two. We do not know the reason for 
this difference. We show how our units convert to those 
in Kochanek & White (2000) as follows: p£ w = 27rp , 
t*J = 1.7 x 2\/2at rfi = llt rfi and ag^f = 2na . Their 
simulation for a^Yi = 1, which has the same cross section 
as ours, reached minimum density at t w — lli r ,n 

and collapsed gravothermally to p c — 2p™ = 13 po at 
t ~ 3.2iJ5* = 35t r ,o, while our simulation reached those 
densities at 20t r ,o and 70i r ,n, respectively (Fig. 6). Their 
evolution is about twice as fast as in our simulation. The 
origin of this discrepancy is unknown. 

4-1-4 The Source of Difference 

Our Monte Carlo JV-body simulations agree with the so- 
lutions of the conducting fluid model in the self-similar 
gravothermal collapse phase, but in general have about a 
20 per cent difference in the central density and the collapse 
rate. This modest disagreement is not surprising in view of 
the fact that the conducting fluid model is not an exact the- 
ory derived from first principals. Heggie & Stevenson (1988) 
compared the thermal conductivity of the conducting fluid 
model with that calculated from the orbit-averaged Fokker- 
Planck equation for star clusters with several profiles, in- 
cluding poly tropes and lowered Maxwellians, evaluated at 
the centres. They find that the coefficients of conductivity 
C varied from one profile to another by factors of 2 or 3. 
The overall collapse rates are not that different, probably 
because the profiles quickly converge to the self-similar so- 
lution around the centre. Indeed, for star clusters, the value 
of C = 0.88, which makes the fluid model match the asymp- 
totic collapse rate in the isotropic Fokker-Planck calcula- 
tion (Cohn 1980), also gives the correct collapse time of the 
Plummer model, 15.4 half-mass relaxation times (Goodman 
1987; Heggie & Ramamani 1989). This means that C needs 
not to be different for the cases of self-similar collapse solu- 
tion and the collapse of Plummer profile. 

One possible distinction between the ID conducting 
fluid model and the 3D iV-body/Monte Carlo simulations 
is that the former assumes that particle velocity distribu- 



tions are isotropic in the frame of the bulk flow while the 
latter do not. Anisotropy in velocity dispersion affects the 
collapse time (Bettwieser 1983; Louis 1990). Fokker-Planck 
calculations with anisotropy show that the collapse time 
of the star clusters initialized by the Plummer model is 
20 per cent larger than for the isotropic case (Takahashi 
1995), and agrees with TV-body simulation (Khalisi et al. 
2007). In Fig. 7, we plot the radial (v%) and tangential 
[v± = (vj + i$)/2] velocity dispersions of our simulation 
when the central density has increased by a factor of 10. 
We do not see anisotropies near the centre, but we do at 
r > 5r c for runs P, NFW and H. When particles scatter 
from the centre to large radii, their nearly radial orbits bring 
anisotropy to the initially isotropic velocity distribution at 
large radii. If the original unscattered particles at large radii 
are less numerous, because the density profile is steeper, the 
anisotropy that results from scattering particles from small 
to large radii will be relatively larger. This is why anisotropy 
is larger if the logarithmic slope of the density profile is 
steeper. Anisotropy may therefore play some role in the col- 
lapse rate of run P and H. 

In short, the conducting fluid model has been shown 
to describe the gravothermal collapse relatively accurately, 
but we cannot expect very high precision in general. There- 
fore, the 20 per cent match in minimum core density of runs 
NFW and H, and 20 per cent match in the gravothermal col- 
lapse time for runs P, NFW and H run is a very reasonable 
agreement. 

4.2 Transitional Regime 

When Kn becomes comparable to, or smaller than, one, the 
self-similar collapse solution of BSI no longer applies, be- 
cause of the presence of the second term in the heat con- 
duction equation (3). The time evolution in units of the 
relaxation time becomes slower the smaller Kn gets, be- 
cause smaller mean free path suppresses particle transport 
more, which results in smaller heat conduction (BSI). To 
test the heat conduction equation (3) when both long and 
short mean-free-path terms are important (i.e. transitional 
regime) , we compared the time evolution of Monte Carlo N- 
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CT C (0) 


TV 


A/r 32 


tio/*r,c(0) 


0.25 


4 x 64 3 


10.3 


374 


0.5 


4 x 64 3 


5.2 


417 


0.75 


2 x 128 3 


5.5 


510 


1.0 


2 x 128 3 


4.1 


585 



10 
8 
6 



Table 2. The transitional regime simulations. a c is the initial 
dimensionless cross section, TV is the number of TV-body particles, 
-V r 32(0) is the ratio of the mean free path to the kernel size at 
the centre at t = 0, and iio is the time for the density to increase 
by a factor of 10. Other parameters are the same as for the BSI 
run. 



body simulations with those of the conducting fluid model, 
derived numerically when no self-similar solution exists, with 
initial Knudsen numbers tfn _1 (0) = <r c (0) = 0.25,0.5,0.75 
and 1.0. The initial condition is the same as in Section 4.1.1: 
the particle distribution of the BSI self-similar collapse solu- 
tion. See Table 2 for the summary. TV-body simulation with 
larger a c becomes very difficult because of the mean free 
path requirement equation (16). The initial ratio of mean 
free path to the kernel size, A/V32 oc p~ 2 / 3 , at the centre is 
listed in the table. Our cr c (0) = 1.0 run violates A > r 3 2 at 
pc(t) ~ 8. The number of particles required to follow the 
collapse scales as beyond this density or cross section. 

For the fluid model, the prefactor C is chosen to be 0.75 
to make an agreement at small a c (Section 4.1.1). We ran the 
fluid code with initial dimensionless cross section a(0)/Vb = 
0,0.25,0.5,0.75, 1.0, 1.5 and 2.0. The time evolution of the 
fluid model in units of initial relaxation time depends only 
on the combination a/vb, which can be seen from the heat 
conduction equation (3). 

We plot the evolution of the central density for both 
TV-body and fluid model results in Fig. 8. In units of the 
relaxation time, the collapse is slower for larger cross sec- 
tion. Fluid model results for a c (Q)/Vb = 0.5, 1.0, 1.5 and 2.0 
evolve similarly to TV-body simulation results for <r c (0) = 
0.25, 0.5, 0.75 and 1.0, respectively. This suggest that the ef- 
fective value of the coefficient is b — 0.25 in the transitional 
regime, Kn > 1. 

To show the agreement between our TV-body results and 
the fluid model solution with b = 0.25, we plotted the nor- 
malized collapse rate as a function of cross section in Fig. 9, 
where fio is the time for the density to increase by a factor 
of 10, and the fiducial time scale t\ is the time at which 
the self-similar solution with <r c (0) = 1 has a density in- 
crease by a factor of 10. The collapse rate is proportional 
to the cross section in the large Kn regime (small a), but 
deviates from the linear relation in the transitional regime 
a > 0.5, as predicted by BSI. The collapse rate should reach 
some maximum at some cross section and then decrease as 
t^Q oc 1 as <j c — > 00. Furthermore, since the value of b can 
be calculated from first principals in the small Kn regime, 
b should converge to the Chapman-Enskog value (b=1.38 in 
Section 2) as Kn — > 0. However, due to the numerical limit 
in equation (16), we cannot go into the small Kn regime to 
see the convergence to the Chapman-Enskog theory or the 
turn over of the collapse rate. Our TV-body results are con- 
sistent with a constant b = 0.25 in the range we are able to 
simulate. 
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Figure 8. Gravothermal collapse in the transitional regime: 
Central density against time from Monte Carlo TV-body simu- 
lations and the conducting fluid model. Cross sections of TV- 
body runs are ct c (0) = 0.25,0.50,0.75 and 1.0, from left to 
right, shifted to the right by 300 for readability. Four copies of 
smooth curves are the solutions of the fluid model with cross sec- 
tions & c (0)/Vb = 0,0.25,0.5,1.0,1.5 and 2.0, from left to right. 
First two curves, <j c /Vb = and 0.25, are indistinguishable, and 
<?c/Vb = 0.5, 1.0,0.15, and 2.0 overlap with the TV-body results 
for <t c = 0.25, 0.5, 0.75, and 1.0, respectively. This match suggests 
b = 0.25. 
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Figure 9. Collapse rate plotted as a function of cross section 
for self-similar solution (dotted line), fluid model with previously 
assumed value b = 1.0 (dashed line), with b = 0.25 (solid line) 
and TV-body (crosses). See text for the definition of the collapse 
rate. A value of b = 0.25 makes the fluid model solution (with its 
conductivity that maps smoothly between the limits of long and 
short mean free path) agree with TV-body simulations even in the 
transitional regime. 



5 DISCUSSION: IMPLICATIONS FOR THE 
COSMOLOGICAL SIMILARITY SOLUTION 

In this section, we discuss the consequences for the cosmo- 
logical similarity solution derived by Ahn & Shapiro (2005, 
hereafter, A&S) of replacing the values of the prefactors (in 
the heat conduction equation 3) C = b = 1.0 by C = 0.75 
and b — 0.25, as calibrated by our Monte Carlo TV-body 
simulations. Until now, we have limited our attention to 
the gravothermal relaxation of isolated SIDM haloes. How- 
ever, the haloes that form in a cosmological context are not 
isolated but rather build up over time from the nonlinear 
growth of small-amplitude initial density perturbations in 
the expanding cosmological background universe, which re- 
sults in a continuous infall of additional mass. To model 
this cosmological formation and evolution of SIDM haloes, 
A&S added the BSI heat conduction term to the fully time- 
dependent conservation equations of the fluid approxima- 
tion in ID, spherical symmetry, as described in A&S. They 
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used these equations to derive an analytical solution for the 
so-called 'secondary infall' problem, in an Einstein-de Sit- 
ter universe, with initial perturbation given by the spherical 
overdensity profile SM/M oc M -1 '' 6 , where M(r) is the un- 
perturbed mass in a sphere of radius r at the mean density 
of the universe, and SM(r) — M(r) — M(r) is the mass 
perturbation inside the radius. The solution with this ini- 
tial overdensity is self-similar, that is, the solution is time- 
independent if radius and density are measured in units 
of the time-varying turn-around radius r ta and background 
critical density pb, respectively. The family of similarity so- 
lutions is parametrized by the dimensionless cross section, 

Q = Pb<rr v ir, (22) 

where r vir is the halo virial radius, the radius at which an 
accretion shock occurs. The density profile has a core whose 
density and size depend upon the value of Q. 

The heat conduction equation (3) has a dependence on 
a = H/X given by, 

L(xa (a- 1 ^ 1 +b- 1 a 2 y 1 . (23) 

This function takes its maximum value y/abC /2 at a = 
Va~ 1 bC~ 1 . Compared to previously assumed values, C = 
b — 1.0, our calibrated heat conduction has a maximum that 
is smaller by a factor of \fb~C « 0.4, which occurs at a value 
of a that is smaller by a factor of V&C -1 = 0.58. There is a 
minimum central density and largest core size which occur 
approximately when the dimensionless cross section at the 
centre, a c , defined in Section 3.2, maximizes the heat flux 
with the value, o c = Va~ 1 bC~ 1 . The solution with that a c 
is defined as the maximally-relaxed solution, and the cor- 
responding cross section Q is denoted by Qth (A&S). This 
maximally-relaxed halo has a density profile almost identi- 
cal to the empirical Burkert profile (Burkert 1995), which 
fits the observed rotation curves of dwarf and LSB galaxies 
well. 

We used the same numerical code as A&S to solve 
for the cosmological similarity solutions with our calibrated 
prefactors, and to compare with their original solutions. 
In Fig. 10, we plot the maximally-relaxed density profile 
(left panel), and the dependence of the central density on 
dimensionless cross section a c (right panel). For prefac- 
tors C = b = 1.0 adopted by A&S, the solution is maxi- 
mally relaxed for Qth = 7.35 x 10~ 4 , with central density 
p c = 1.17xl0 4 p 6 . ForC = 0.75 and b = 0.25, the maximally- 
relaxed solution shifts to Qth = 3.95 x 10 -4 , with a central 
density p c = 1.50 x 10 4 pb. 5 This central density is about 
30 per cent larger than that of the original solution. For 
a c <S 1, a cross section that is 1.33 times as large is required 
to achieve the same central density, due to the change of 
coefficient C from 1.0 to 0.75. 



5 The value of u c which defined the maximally-relaxed solution, 
&c = Va~ 1 bC~ 1 , changes from 0.666 to 0.384 if the conductivity 
parameters change their values in A&S to our new values here. 
However, the actual minimum densities occur at slightly higher 
a. For C = b = 1.0, central density takes its minimum value 
p c = 1.15 x 10 4 p 6 at a c = 0.85, or Q = 9.78 x 10~ 4 . For C = 0.75 
and b = 0.25, the minimum density p c = 1.46 X 10 4 p(, is found at 
ct c = 0.55, or Q = 3.95 x 10~ 4 . 




Figure 10. The cosmological similarity solution by Ahn & 
Shapiro (2005) with our calibrated prefactors (C=0.75, b=0.25; 
solid line) and the original prefactors (C=b=1.0; dotted line). 
(Left:) The maximally-relaxed density profiles, with density and 
radius in units of background critical density pb and turn around 
radius rta- The maximally-relaxed solution for new prefactors has 
a 30 per cent larger central density. (Right:) The central density 
for old and new prefactors for different cross section values. 



A&S estimated the cross section values that brought 
dwarf galaxy rotation curves into agreement with the cos- 
mological self-similar halo profiles. In order to translate a 
given value of Qth into a cross section value a, the typical 
formation time for haloes that host dwarf galaxies in the 
ACDM model of structure formation was used to relate r v ir 
and pb in the definition of Q, as described in A&S. Using 
the same argument, our corrected Qth value corresponds to 
a cross section 117 cm 2 g" 1 , which is smaller than the origi- 
nal value 218 cm 2 g _1 but still much larger than the values, 
0.5 — 5 cm 2 g _1 , found in Monte Carlo iV-body simulations 
that produce cored profiles for galactic haloes formed from 
cosmological initial conditions (Dave et al. 2001). 

We will discuss this apparent discrepancy between the 
values of the SIDM cross section which are best able to 
match the observed rotation curves of dwarf and LSB galax- 
ies, from the A&S similarity solutions and cosmological N- 
body/Monte Carlo simulations, respectively, in a separate 
paper. Here, we have demonstrated that this discrepancy is 
not likely to result from the break down of the conducting 
fluid model. We have found good agreement, in fact, between 
the Monte Carlo A^-body simulations and the conducting 
fluid model for the thermal relaxation and gravothermal 
collapse of isolated haloes (of fixed mass), at least with re- 
gard to the long mean-free-path regime and the transitional 
regime in which the mean free path is comparable to the 
system size. These are the regimes of greatest relevance to 
the SIDM halo problem. We have also shown that, when we 
use the agreement between our Monte Carlo A^-body results 
and the conducting fluid model solutions for isolated haloes 
to calculate the dimensionless parameters on which the heat 
conduction depends, the impact of the modified parameters 
on the A&S similarity solution is relatively small. We must, 
therefore, seek a different explanation for the different val- 
ues of a required by the cosmological A"-body/Monte Carlo 
and the similarity solutions of A&S, to produce the same 
observed degree of relaxation of the SIDM halo density pro- 
files. 

As we shall discuss in our next paper in this series, the 
essential difference between the self-similar solution of A&S 
and the cosmological A-body simulations is that the self- 
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similar halo is continuously heated by the accretion shock, 
while the inner part of a more realistic halo, on the galaxy 
scale, is eventually unaffected by infall when infall rate slows 
after the halo forms. We will compare cosmological Monte 
Carlo iV-body simulations with the conducting fluid model 
with non-self-similar infall in our next paper. 



panion paper, in fact, the discrepancy results, instead, from 
the gradual departure of halo evolution from self-similarity 
as the infall rate during cosmological structure formation 
drops below the self-similar rate at late times after halo for- 
mation. 



6 CONCLUSION AND SUMMARY 

The ability of the SIDM hypothesis to resolve the cusp-core 
problem of collisionless CDM haloes on the galaxy scale and 
the cross section values required to accomplish this remain 
uncertain, as long as the previous conclusions drawn from 
./V-body simulation and the conducting fluid model differ 
so strongly. To rectify this situation, we have developed a 
new Monte Carlo iV-body code of our own, based on the 
pre-existing GADGET 1.1 iV-body code, and applied it to 
compare the iV-body simulations with the conducting fluid 
model for isolated, spherically symmetric self-gravitating 
SIDM haloes. The collisions were assumed to be velocity 
independent, elastic, and isotropic. Our results include the 
following: 

(i) Our Monte Carlo iV-body simulations are in very good 
agreement with the analytical self-similar gravothermal col- 
lapse solution of BSI, when the coefficient of thermal con- 
duction is set to C — 0.75 (Section 4.1.1); the density and 
velocity dispersion profiles evolve self-similar ly; the central 
density and velocity dispersion follow the self-similar time 
evolution formulae with the predicted constant a = 2.19; 
and, the time to collapse is always proportional to the cen- 
tral relaxation time at that time (equation 6-10). 

(ii) The conducting fluid model agrees with Monte Carlo 
./V-body simulations reasonably well for different initial con- 
ditions: Plummer model, Hernquist profile, and NFW pro- 
file, in the large- Kn regime. The collapse time and the cen- 
tral density at maximum core expansion agree within 20 per 
cent. The shape of the density profile and the central density 
evolution as a function of time during gravothermal collapse 
agree very well. 

(iii) We also showed that the collapse time becomes longer 
in units of relaxation time as the system transitions from the 
large- to the small- Kn regime, as predicted by the conduct- 
ing fluid model. The iV-body results agree with the con- 
ducting fluid model for Kn 1, or a ^ 1, with the prefactor 
b — 0.25. However, this prefactor is more than five times 
smaller than the Chapman-Enskog value, valid asymptoti- 
cally in the small Kn limit. The conducting fluid model must 
be further calibrated against iV-body simulations if it is used 
beyond the transitional regime, for a > 1. 

(iv) Our calibration of the prefactors C and b does 
not change the cosmological similarity solutions of Ahn & 
Shapiro (2005) significantly. The cross section that gives the 
minimum central density on the dwarf galaxy scale is al- 
tered from 220 cm 2 g _1 to 117 cm 2 g _1 , but this is still much 
larger than the values that cosmological Monte Carlo N- 
body simulations used to make cored SIDM haloes (a ~ 
0.5 — 5 cm 2 g _1 ). We will investigate this problem in our 
subsequent paper. Our results here suggest that this appar- 
ent discrepancy is not the result of a breakdown of either the 
conducting fluid model or the Monte Carlo scattering algo- 
rithm in the JV-body simulations. As we shall show in a com- 
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